home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / pairalign.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  14.1 KB  |  584 lines  |  [TEXT/R*ch]
  1. /* Change int h to int gh everywhere  DES June 1994 */
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <string.h>
  5. #include <math.h>
  6. #include "clustalw.h"
  7.  
  8. #define MIN(a,b) ((a)<(b)?(a):(b))
  9. #define MAX(a,b) ((a)>(b)?(a):(b))
  10.  
  11. #define gap(k)  ((k) <= 0 ? 0 : g + gh * (k))
  12. #define tbgap(k)  ((k) <= 0 ? 0 : tb + gh * (k))
  13. #define tegap(k)  ((k) <= 0 ? 0 : te + gh * (k))
  14.  
  15.  
  16. /*
  17. *    Prototypes
  18. */
  19.  
  20. extern void *ckalloc(size_t);
  21. extern void *ckfree(void *);
  22. extern int  get_matrix(int *matptr, int *xref, int matrix[NUMRES][NUMRES], int pos_flag);
  23.  
  24. void pairalign(int istart, int iend, int jstart, int jend);
  25. static double     tracepath(short sb1, short sb2);
  26. static void     add(short);
  27. static void     del(short);
  28. static int       calc_score(short,short,short,short);
  29. static int       diff(short,short,short,short,short,short);
  30. static void     forward_pass(char *ia, char *ib, short n, short m);
  31. static void     reverse_pass(char *ia, char *ib, short n, short m);
  32.  
  33. /*
  34.  *   Global variables
  35.  */
  36.  
  37. extern double   **tmat;
  38. extern float    pw_go_penalty;
  39. extern float    pw_ge_penalty;
  40. extern int     nseqs;
  41. extern int     max_aa;
  42. extern int     *seqlen_array;
  43. extern int     debug;
  44. extern int      mat_avscore;
  45. extern int     blosum30mt[], pam350mt[], idmat[], pw_usermat[];
  46. extern int    blosum45mt[], md_350mt[];
  47. extern int    def_aa_xref[], pw_aa_xref[];
  48. extern Boolean  dnaflag, is_weight;
  49. extern char     **seq_array;
  50. extern char     *amino_acid_codes;
  51. extern char     pw_mtrxname[];
  52.  
  53. static double     mm_score;
  54. static int     print_ptr,last_print;
  55. /* static int     displ[2*MAXLEN];   */
  56. /* static int     HH[MAXLEN+1], DD[MAXLEN+1], RR[MAXLEN+1], SS[MAXLEN+1]; */
  57. extern int *displ, *HH, *DD, *RR, *SS;   /* dgg made these extern */
  58.  
  59. static short     g, gh;
  60. static int       seq1, seq2;
  61. static int        matrix[NUMRES][NUMRES];
  62. static int        maxscore;
  63. static short        sb1, sb2, se1, se2;
  64.  
  65. void pairalign(int istart, int iend, int jstart, int jend)
  66. {
  67.   int     *mat_xref;
  68.   static short    si, sj, i, j, k;
  69.   static short    n,m,len1,len2;
  70.   static short    length, itemp, end = FALSE;
  71.   static short    maxres;
  72.   static int    *matptr;
  73.   static char   c;
  74.   static int    sum;
  75.   char matfile[FILENAMELEN];
  76.  
  77.   if (dnaflag)
  78.      {
  79.        for(i=0;i<5;++i)
  80.          {
  81.             for(j=0;j<5;++j) {
  82.                  if(i==j)
  83.                      matrix[i][j]=1000;
  84.                  else
  85.                      matrix[j][i]=0;
  86.              }
  87.             if(is_weight) {
  88.                 matrix[0][2]=500;
  89.                 matrix[2][0]=500;
  90.                 matrix[1][3]=500;
  91.                 matrix[3][1]=500;
  92.                 matrix[1][4]=500;
  93.                 matrix[4][1]=500;
  94.              }
  95.          }
  96.     }
  97.   else
  98.     {
  99. if (debug>1) fprintf(stderr,"matrix %s\n",pw_mtrxname);
  100.        if (strcmp(pw_mtrxname, "blosum") == 0)
  101.           {
  102.              matptr = blosum30mt;
  103.              mat_xref = def_aa_xref;
  104.           }
  105.        else if (strcmp(pw_mtrxname, "md") == 0)
  106.           {
  107.              matptr = md_350mt;
  108.              mat_xref = def_aa_xref;
  109.           }
  110.        else if (strcmp(pw_mtrxname, "pam") == 0)
  111.           {
  112.              matptr = pam350mt;
  113.              mat_xref = def_aa_xref;
  114.           }
  115.        else if (strcmp(pw_mtrxname, "id") == 0)
  116.           {
  117.              matptr = idmat;
  118.              mat_xref = def_aa_xref;
  119.           }
  120.        else
  121.           {
  122.              matptr = pw_usermat;
  123.              mat_xref = pw_aa_xref;
  124.           }
  125.  
  126.        maxres = get_matrix(matptr, mat_xref, matrix, TRUE);
  127.        if (maxres == 0) exit(-1);
  128.     }
  129.  
  130.  
  131.   for (si=MAX(0,istart);si<nseqs && si<iend;si++)
  132.    {
  133.      n = seqlen_array[si+1];
  134.      len1 = 0;
  135.      for (i=1;i<=n;i++) {
  136.     c = seq_array[si+1][i];
  137.     if ((c>=0) && (c <= max_aa)) len1++;
  138.      }
  139.  
  140.      for (sj=MAX(si+1,jstart+1);sj<nseqs && sj<jend;sj++)
  141.       {
  142.         m = seqlen_array[sj+1];
  143.     len2 = 0;
  144.     for (i=1;i<=m;i++) {
  145.         c = seq_array[sj+1][i];
  146.         if ((c>=0) && (c <= max_aa)) len2++;
  147.     }
  148.  
  149. if (n>MAXLEN) fprintf(stderr,"Error: MAXLEN too small, seq 1 is %d long\n",(pint)n);
  150. if (m>MAXLEN) fprintf(stderr,"Error: MAXLEN too small, seq 1 is %d long\n",(pint)m);
  151.  
  152.         if (dnaflag) {
  153.            g = 200.0 * (float)pw_go_penalty;
  154.            gh = pw_ge_penalty * 100.0;
  155.         }
  156.         else {
  157.            if (mat_avscore <= 0)
  158.               g = 200.0 * (float)(pw_go_penalty + log((double)(MIN(n,m))));
  159.            else
  160.               g = 2.0 * mat_avscore * (float)(pw_go_penalty +
  161.                     log((double)(MIN(n,m))));
  162.            gh = pw_ge_penalty * 100.0;
  163.         }
  164.  
  165. if (debug>1) fprintf(stderr,"go %d ge %d\n",(pint)g,(pint)gh);
  166.  
  167. /*
  168.    align the sequences
  169. */
  170.         seq1 = si+1;
  171.         seq2 = sj+1;
  172.  
  173.         forward_pass(&seq_array[seq1][0], &seq_array[seq2][0],
  174.            n, m);
  175.  
  176.         reverse_pass(&seq_array[seq1][0], &seq_array[seq2][0],
  177.            n, m);
  178.  
  179.         last_print = 0;
  180.     print_ptr = 1;
  181. /*
  182.         sb1 = sb2 = 1;
  183.         se1 = n-1;
  184.         se2 = m-1;
  185. */
  186.  
  187. /* use Myers and Miller to align two sequences */
  188.  
  189.         maxscore = diff(sb1-1, sb2-1, se1-sb1+1, se2-sb2+1, 
  190.         (short)0, (short)0);
  191.  
  192. /* calculate percentage residue identity */
  193.  
  194.         mm_score = tracepath(sb1,sb2);
  195.  
  196.     mm_score /= (double)MIN(len1,len2);
  197.  
  198.         tmat[si+1][sj+1] = ((double)100.0 - mm_score)/(double)100.0;
  199.         tmat[sj+1][si+1] = ((double)100.0 - mm_score)/(double)100.0;
  200.  
  201. if (debug>0)
  202. {
  203.         fprintf(stdout,"Sequences (%d:%d) Aligned. Score: %d CompScore:  %d\n",
  204.                            (pint)si+1,(pint)sj+1, 
  205.                            (pint)mm_score, 
  206.                            (pint)maxscore/(MIN(len1,len2)*100));
  207. }
  208. else
  209. {
  210.         fprintf(stdout,"Sequences (%d:%d) Aligned. Score:  %d\n",
  211.                                       (pint)si+1,(pint)sj+1, 
  212.                                       (pint)mm_score);
  213. }
  214.  
  215.    }
  216.   }
  217.  
  218.   return;
  219. }
  220.  
  221. static void add(short v)
  222. {
  223.  
  224.         if(last_print<0) {
  225.                 displ[print_ptr-1] = v;
  226.                 displ[print_ptr++] = last_print;
  227.         }
  228.         else
  229.                 last_print = displ[print_ptr++] = v;
  230. }
  231.  
  232. static int calc_score(short iat,short jat,short v1,short v2)
  233. {
  234.         short ipos,jpos;
  235.     int ret;
  236.  
  237.         ipos = v1 + iat;
  238.         jpos = v2 + jat;
  239.  
  240.         ret=matrix[seq_array[seq1][ipos]][seq_array[seq2][jpos]];
  241.  
  242.     return(ret);
  243. }
  244.  
  245.  
  246. static double tracepath(short tsb1, short tsb2)
  247. {
  248.     char c1,c2;
  249.         short  i1,i2;
  250.         short i,j,k,r,pos,to_do;
  251.     short count, total;
  252.     double score;
  253. char s1[MAXLEN*2+1], s2[MAXLEN*2+1];
  254.  
  255. if ((print_ptr>MAXLEN*2)||(print_ptr<0))
  256.  fprintf(stderr,"Error: displ array too small %d\n",(pint)print_ptr);
  257.         to_do=print_ptr-1;
  258.         i1 = tsb1;
  259.         i2 = tsb2;
  260.  
  261.     pos = 0;
  262.     count = 0;
  263.         for(i=1;i<=to_do;++i) {
  264. if ((i1>MAXLEN)||(i1<0)) fprintf(stderr,"Error: index i1 %d\n",(pint)i1);
  265. if ((i2>MAXLEN)||(i2<0)) fprintf(stderr,"Error: index i2 %d\n",(pint)i2);
  266. if ((pos>MAXLEN*2)||(pos<0)) fprintf(stderr,"Error: index pos %d\n",(pint)pos);
  267.  
  268. if (debug>1) fprintf(stderr,"%d ",(pint)displ[i]);
  269.                 if(displ[i]==0) {
  270.             c1 = seq_array[seq1][i1];
  271.             c2 = seq_array[seq2][i2];
  272. if (debug>1)
  273. {
  274. if (c1>max_aa) s1[pos] = '-';
  275. else s1[pos]=amino_acid_codes[c1];
  276. if (c2>max_aa) s2[pos] = '-';
  277. else s2[pos]=amino_acid_codes[c2];
  278. }
  279.             if ((c1>=0) && (c1 <= max_aa) &&
  280.                                     (c1 == c2)) count++;
  281.                         ++i1;
  282.                         ++i2;
  283.                         ++pos;
  284.                 }
  285.                 else {
  286.                         if((k=displ[i])>0) {
  287. if (debug>1)
  288. for (r=0;r<k;r++)
  289. {
  290. s1[pos+r]='-';
  291. if (seq_array[seq2][i2+r]>max_aa) s2[pos+r] = '-';
  292. else s2[pos+r]=amino_acid_codes[seq_array[seq2][i2+r]];
  293. }
  294.                                 i2 += k;
  295.                                 pos += k;
  296.                                 count--;
  297.                         }
  298.                         else {
  299. if (debug>1)
  300. for (r=0;r<(-k);r++)
  301. {
  302. s2[pos+r]='-';
  303. if (seq_array[seq1][i1+r]>max_aa) s1[pos+r] = '-';
  304. else s1[pos+r]=amino_acid_codes[seq_array[seq1][i1+r]];
  305. }
  306.                                 i1 -= k;
  307.                                 pos -= k;
  308.                                 count--;
  309.                         }
  310.                 }
  311.         }
  312. if (debug>1) fprintf(stderr,"\n");
  313. if (debug>1) 
  314. {
  315. for (i=0;i<pos;i++) fprintf(stderr,"%c",s1[i]);
  316. fprintf(stderr,"\n");
  317. for (i=0;i<pos;i++) fprintf(stderr,"%c",s2[i]);
  318. fprintf(stderr,"\n");
  319. }
  320.  
  321.         if (count <= 0) count = 1;
  322.     score = (double)100.0 * (double)count;
  323.     return(score);
  324. }
  325.  
  326.  
  327. static void forward_pass(char *ia, char *ib, short n, short m)
  328. {
  329.  
  330.   short i,j,k;
  331.   int f,hh,p,t;
  332.  
  333.   maxscore = 0;
  334.   se1 = se2 = 0;
  335.   for (i=0;i<=m;i++)
  336.     {
  337.        HH[i] = 0;
  338.        DD[i] = -g;
  339.     }
  340.  
  341.   for (i=1;i<=n;i++)
  342.      {
  343.         hh = p = 0;
  344.     f = -g;
  345.  
  346.         for (j=1;j<=m;j++)
  347.            {
  348.  
  349.               f -= gh; 
  350.               t = hh - g - gh;
  351.               if (f<t) f = t;
  352.  
  353.               DD[j] -= gh;
  354.               t = HH[j] - g - gh;
  355.               if (DD[j]<t) DD[j] = t;
  356.  
  357.               hh = p + matrix[ia[i]][ib[j]];
  358.               if (hh<f) hh = f;
  359.               if (hh<DD[j]) hh = DD[j];
  360.               if (hh<0) hh = 0;
  361.  
  362.               p = HH[j];
  363.               HH[j] = hh;
  364.  
  365.               if (hh > maxscore)
  366.                 {
  367.                    maxscore = hh;
  368.                    se1 = i;
  369.                    se2 = j;
  370.                 }
  371.            }
  372.      }
  373.  
  374. }
  375.  
  376.  
  377. static void reverse_pass(char *ia, char *ib, short n, short m)
  378. {
  379.  
  380.   short i,j,k;
  381.   int f,hh,p,t;
  382.   int cost;
  383.  
  384.   cost = 0;
  385.   sb1 = sb2 = 0;
  386.   for (i=se2;i>0;i--)
  387.     {
  388.        HH[i] = -1;
  389.        DD[i] = -1;
  390.     }
  391.  
  392.   for (i=se1;i>0;i--)
  393.      {
  394.         hh = f = -1;
  395.         if (i == se1) p = 0;
  396.         else p = -1;
  397.  
  398.         for (j=se2;j>0;j--)
  399.            {
  400.  
  401.               f -= gh; 
  402.               t = hh - g - gh;
  403.               if (f<t) f = t;
  404.  
  405.               DD[j] -= gh;
  406.               t = HH[j] - g - gh;
  407.               if (DD[j]<t) DD[j] = t;
  408.  
  409.               hh = p + matrix[ia[i]][ib[j]];
  410.               if (hh<f) hh = f;
  411.               if (hh<DD[j]) hh = DD[j];
  412.  
  413.               p = HH[j];
  414.               HH[j] = hh;
  415.  
  416.               if (hh > cost)
  417.                 {
  418.                    cost = hh;
  419.                    sb1 = i;
  420.                    sb2 = j;
  421.                    if (cost >= maxscore) break;
  422.                 }
  423.            }
  424.         if (cost >= maxscore) break;
  425.      }
  426.  
  427. }
  428.  
  429. static int diff(short A,short B,short M,short N,short tb,short te)
  430. {
  431.         short midi,midj,i,j,type;
  432.         int midh;
  433.         static int f, hh, c, e, s, t;
  434.  
  435.         if(N<=0)  {
  436.                 if(M>0) {
  437.                         del(M);
  438.                 }
  439.  
  440.                 return(-tbgap(M));
  441.         }
  442.  
  443.         if(M<=1) {
  444.                 if(M<=0) {
  445.                         add(N);
  446.                         return(-tbgap(N));
  447.                 }
  448.  
  449.                 midh = -(tb+gh) - tegap(N);
  450.                 hh = -(te+gh) - tbgap(N);
  451.         if (hh>midh) midh = hh;
  452.                 midj = 0;
  453.                 for(j=1;j<=N;j++) {
  454.                         hh = calc_score(1,j,A,B)
  455.                                     - tegap(N-j) - tbgap(j-1);
  456.                         if(hh>midh) {
  457.                                 midh = hh;
  458.                                 midj = j;
  459.                         }
  460.                 }
  461.  
  462.                 if(midj==0) {
  463.                         del(1);
  464.                         add(N);
  465.                 }
  466.                 else {
  467.                         if(midj>1)
  468.                                 add(midj-1);
  469.                         displ[print_ptr++] = last_print = 0;
  470.                         if(midj<N)
  471.                                 add(N-midj);
  472.                 }
  473.                 return midh;
  474.         }
  475.  
  476. /* Divide: Find optimum midpoint (midi,midj) of cost midh */
  477.  
  478.         midi = M / 2;
  479. if (N>MAXLEN) fprintf(stderr,"Error: MAXLEN too small, seq 1 is %d long\n",(pint)N);
  480. if (M>MAXLEN) fprintf(stderr,"Error: MAXLEN too small, seq 2 is %d long\n",(pint)M);
  481.         HH[0] = 0.0;
  482.         t = -tb;
  483.         for(j=1;j<=N;j++) {
  484.                 HH[j] = t = t-gh;
  485.                 DD[j] = t-g;
  486.         }
  487.  
  488.         t = -tb;
  489.         for(i=1;i<=midi;i++) {
  490.                 s=HH[0];
  491.                 HH[0] = hh = t = t-gh;
  492.                 f = t-g;
  493.                 for(j=1;j<=N;j++) {
  494.                         if ((hh=hh-g-gh) > (f=f-gh)) f=hh;
  495.                         if ((hh=HH[j]-g-gh) > (e=DD[j]-gh)) e=hh;
  496.                         hh = s + calc_score(i,j,A,B);
  497.                         if (f>hh) hh = f;
  498.                         if (e>hh) hh = e;
  499.  
  500.                         s = HH[j];
  501.                         HH[j] = hh;
  502.                         DD[j] = e;
  503.                 }
  504.         }
  505.  
  506.         DD[0]=HH[0];
  507.  
  508.         RR[N]=0;
  509.         t = -te;
  510.         for(j=N-1;j>=0;j--) {
  511.                 RR[j] = t = t-gh;
  512.                 SS[j] = t-g;
  513.         }
  514.  
  515.         t = -te;
  516.         for(i=M-1;i>=midi;i--) {
  517.                 s = RR[N];
  518.                 RR[N] = hh = t = t-gh;
  519.                 f = t-g;
  520.  
  521.                 for(j=N-1;j>=0;j--) {
  522.  
  523.                         if ((hh=hh-g-gh) > (f=f-gh)) f=hh;
  524.                         if ((hh=RR[j]-g-gh) > (e=SS[j]-gh)) e=hh;
  525.                         hh = s + calc_score(i+1,j+1,A,B);
  526.                         if (f>hh) hh = f;
  527.                         if (e>hh) hh = e;
  528.  
  529.                         s = RR[j];
  530.                         RR[j] = hh;
  531.                         SS[j] = e;
  532.  
  533.                 }
  534.         }
  535.  
  536.         SS[N]=RR[N];
  537.  
  538.         midh=HH[0]+RR[0];
  539.         midj=0;
  540.         type=1;
  541.         for(j=0;j<=N;j++) {
  542.                 hh = HH[j] + RR[j];
  543.                 if(hh>=midh)
  544.                         if(hh>midh || HH[j]!=DD[j] && RR[j]==SS[j]) {
  545.                                 midh=hh;
  546.                                 midj=j;
  547.                         }
  548.         }
  549.  
  550.         for(j=N;j>=0;j--) {
  551.                 hh = DD[j] + SS[j] + g;
  552.                 if(hh>midh) {
  553.                         midh=hh;
  554.                         midj=j;
  555.                         type=2;
  556.                 }
  557.         }
  558.  
  559.         /* Conquer recursively around midpoint  */
  560.  
  561.  
  562.         if(type==1) {             /* Type 1 gaps  */
  563.                 diff(A,B,midi,midj,tb,g);
  564.                 diff(A+midi,B+midj,M-midi,N-midj,g,te);
  565.         }
  566.         else {
  567.                 diff(A,B,midi-1,midj,tb,0.0);
  568.                 del(2);
  569.                 diff(A+midi+1,B+midj,M-midi-1,N-midj,0.0,te);
  570.         }
  571.  
  572.         return midh;       /* Return the score of the best alignment */
  573. }
  574.  
  575. static void del(short k)
  576. {
  577.         if(last_print<0)
  578.                 last_print = displ[print_ptr-1] -= k;
  579.         else
  580.                 last_print = displ[print_ptr++] = -(k);
  581. }
  582.  
  583.  
  584.